// file: cs-estimate.do
// purpose: replicate main analyses in paper


  quietly log
  local logon = r(status)
  if "`logon'" == "on" { 
  log close 
  }
log using cs-estimate.log, text replace

use cs-recode.dta, clear
 
// Run multilevel models for Tables 2 and 3
 
quietly meologit opst lncs i.year ||fips:

	estimates store m1

quietly meologit opst lncs opfed c.percapincreal unemployment debt_net_pc_real pct_ba ///
	pct_black pct_hisp pct_rural mds1 mds2_abs h_diffs c.policy##c.partyid ///
	r_male c.r_age##c.r_age r_educ r_hisp r_black ///
	c.hou_majority##c.partyid c.unified_gov##c.partyid i.year ||fips:
	
estimates store m2
estimates save m2.est, replace

quietly meologit opst c.lnpop##c.size_low opfed c.percapincreal unemployment debt_net_pc_real pct_ba ///
	pct_black pct_hisp pct_rural mds1 mds2_abs h_diffs  ///
	r_male c.r_age##c.r_age r_educ r_hisp r_black ///
	c.hou_majority##c.partyid c.policy##c.partyid c.unified_gov##c.partyid i.year ||fips:
	
estimates store m4
estimates save m4.est, replace	


quietly meologit opst lncs opfed c.percapincreal unemployment debt_net_pc_real pct_ba ///
	pct_black pct_hisp pct_rural mds1 mds2_abs h_diffs c.policy##c.partyid ///
	r_male c.r_age##c.r_age r_educ r_hisp r_black ///
	c.hou_majority##c.partyid c.unified_gov##c.partyid i.year if stateabv~="NH" & stateabv~="CA" ||fips:
	
estimates store m3
estimates save m3.est, replace

// Now for weighted, single-level models:

ologit opst lncs i.year [pw=weight], cluster(fips)

	estimates store m1_wgt
	
ologit opst lncs opfed c.percapincreal unemployment debt_net_pc_real pct_ba ///
	pct_black pct_hisp pct_rural mds1 mds2_abs h_diffs c.policy##c.partyid ///
	r_male c.r_age##c.r_age r_educ r_hisp r_black ///
	c.hou_majority##c.partyid c.unified_gov##c.partyid i.year [pw=weight], cluster(fips)
	
estimates store m2_wgt	

ologit opst lncs opfed c.percapincreal unemployment debt_net_pc_real pct_ba ///
	pct_black pct_hisp pct_rural mds1 mds2_abs h_diffs c.policy##c.partyid ///
	r_male c.r_age##c.r_age r_educ r_hisp r_black ///
	c.hou_majority##c.partyid c.unified_gov##c.partyid i.year if stateabv~="NH" & stateabv~="CA" [pw=weight], cluster(fips)
	
estimates store m3_wgt

ologit opst c.lnpop##c.size_low opfed c.percapincreal unemployment debt_net_pc_real pct_ba ///
	pct_black pct_hisp pct_rural mds1 mds2_abs h_diffs  ///
	r_male c.r_age##c.r_age r_educ r_hisp r_black ///
	c.hou_majority##c.partyid c.policy##c.partyid c.unified_gov##c.partyid i.year [pw=weight], cluster(fips)
	
estimates store m4_wgt


// Make regression tables using esttab //

label var opfed "Opinions of Federal Government"

esttab m1_wgt m2_wgt m1 m2 using table2_3.txt, cells("b(fmt(3) star) se(fmt(3) par)") ///
	starlevels(+ .1 * .05 ** .01) tex stats(N aic bic) style(tex) ///
	label nolz drop(199*.year 20*.year) ///
	order(lncs policy unified_gov ///
	hou_majority h_diffs mds1 mds2_abs unemployment percapincreal debt_net_pc_real ///
	pct_ba pct_black pct_hisp pct_rural opfed partyid c.policy#c.partyid c.unified_gov#c.partyid ///
	c.hou_majority#c.partyid  ///
	r_male r_educ r_age c.r_age#c.r_age r_black r_hisp)  replace
	
esttab m3_wgt m4_wgt m3 m4 using table2_3.txt, cells("b(fmt(3) star) se(fmt(3) par)") ///
	starlevels(+ .1 * .05 ** .01) tex stats(N aic bic) style(tex) ///
	label nolz drop(199*.year 20*.year) ///
	order(lncs lnpop size_low c.lnpop#c.size_low policy unified_gov ///
	hou_majority h_diffs mds1 mds2_abs unemployment percapincreal debt_net_pc_real ///
	pct_ba pct_black pct_hisp pct_rural opfed partyid c.policy#c.partyid c.unified_gov#c.partyid ///
	c.hou_majority#c.partyid  ///
	r_male r_educ r_age c.r_age#c.r_age r_black r_hisp)  append


// Run robustness models for online appendix //

do cs-alt-measures.do 
do cs-alt-legprof.do	

// Estimate marginal effects and save to dataset //	
	
forvalues z=1/2 {
gen estimate`z' = .
gen se`z' = .

est restore m`z'

margins, expression(predict(pr outcome(3)) + predict(pr outcome(4))) at(lncs = (-1.213 (.2) 3.86)) atmeans post	
margins, coeflegend

replace estimate`z' = _b[1bn._at] in 1
replace se`z' = _se[1bn._at] in 1

forvalues x=2/26 {
replace estimate`z' = _b[`x'._at] in `x'
replace se`z' = _se[`x'._at] in `x'
}
gen lo95_`z' = estimate`z' - 1.96*se`z'
gen up95_`z' = estimate`z' + 1.96*se`z'
}

gen predlncs = .
replace predlncs = -1.213 + .2*(_n-1) if _n<27

bysort fips year: gen numb=_n


// Use predictions to graph marginal effect (Fig 5.)

twoway rarea lo95_2 up95_2 predlncs, sort  fcolor(gs8*.30) lcolor(gs8*.25) || ///
	line estimate2 predlncs, sort yaxis(1) color(black) || ///
	kdensity lncs if numb==1,  color(gs10) lwidth(thin) yaxis(2) ///
	yscale(off axis(2)) scheme(s1mono) plotregion(color(none)) ///
	legend(off) ytitle(Probability of Favorable Opinion ({it:y{subscript:i} {&gt} 2}), ///
	size(small) margin(small)) ///
	xtitle(Constituency Size (log), size(small) margin(small)) ///
	ylabel(,angle(0) labgap(1.25) labsize(small)) ///
	xlabel(, labgap(1.25) labsize(small))
	
graph export marginaleffect.pdf, replace

preserve

keep estimate1-estimate2 se1-se2 lo95* up95* predlncs
drop if predlncs==.

save estimation_results.dta, replace

restore

foreach z in 2_wgt {
gen estimate`z' = .
gen se`z' = .

est restore m`z'

margins, expression(predict(pr outcome(3)) + predict(pr outcome(4))) at(lncs = (-1.213 (.2) 3.86)) atmeans post	
margins, coeflegend

replace estimate`z' = _b[1bn._at] in 1
replace se`z' = _se[1bn._at] in 1

forvalues x=2/26 {
replace estimate`z' = _b[`x'._at] in `x'
replace se`z' = _se[`x'._at] in `x'
}
gen lo95_`z' = estimate`z' - 1.96*se`z'
gen up95_`z' = estimate`z' + 1.96*se`z'
}

gen predlncswgt = .
replace predlncswgt = -1.213 + .2*(_n-1) if _n<27

preserve

keep estimate2_wgt se2_wgt lo95_2_wgt up95_2_wgt predlncswgt
drop if predlncswgt==.

save estimation_results_wgt.dta, replace

restore

foreach z in 2_b {
gen estimate`z' = .
gen se`z' = .

est restore m`z'

margins, expression(predict(pr outcome(3)) + predict(pr outcome(4))) at(lncssen = (.266 (.2) 4.55)) atmeans post	
margins, coeflegend

replace estimate`z' = _b[1bn._at] in 1
replace se`z' = _se[1bn._at] in 1

forvalues x=2/22 {
replace estimate`z' = _b[`x'._at] in `x'
replace se`z' = _se[`x'._at] in `x'
}
gen lo95_`z' = estimate`z' - 1.96*se`z'
gen up95_`z' = estimate`z' + 1.96*se`z'
}

gen predlncssen = .
replace predlncssen = .266 + .2*(_n-1) if _n<23

preserve

keep estimate2_b se2_b lo95_2_b up95_2_b predlncssen
drop if predlncssen==.

save estimation_results_sen.dta, replace

restore


foreach z in 2_c {
gen estimate`z' = .
gen se`z' = .

est restore m`z'

margins, expression(predict(pr outcome(3)) + predict(pr outcome(4))) at(lncsbicam = (-1.27 (.2) 3.45)) atmeans post	
margins, coeflegend

replace estimate`z' = _b[1bn._at] in 1
replace se`z' = _se[1bn._at] in 1

forvalues x=2/24 {
replace estimate`z' = _b[`x'._at] in `x'
replace se`z' = _se[`x'._at] in `x'
}
gen lo95_`z' = estimate`z' - 1.96*se`z'
gen up95_`z' = estimate`z' + 1.96*se`z'
}

gen predlncsbicam = .
replace predlncsbicam = -1.27 + .2*(_n-1) if _n<25

*bysort fips year: gen numb=_n

preserve

keep estimate2_c se2_c lo95_2_c up95_2_c predlncsbicam
drop if predlncsbicam==.

save estimation_results_bicam.dta, replace

restore


// Use results from Table 3, Model 12 and generate estimated change in 
//	predicted effect of increasing chamber size from 75 to 150 seats (Table 4)

*for states of 1 million residents
estimates restore m4
margins, expression(predict(pr outcome(3)) + predict(pr outcome(4)))  at(lnpop=0 size_low=(75 150)) atmeans post pwcompare
mat list e(b)

*for states of 4 million residents
estimates restore m4
margins, expression(predict(pr outcome(3)) + predict(pr outcome(4)))  at(lnpop=1.3862944 size_low=(75 150)) atmeans post pwcompare
mat list e(b)

*for states of 10 million residents
estimates restore m4
margins, expression(predict(pr outcome(3)) + predict(pr outcome(4)))  at(lnpop=2.3025851 size_low=(75 150)) atmeans post pwcompare
mat list e(b)

*for states of 20 million residents
estimates restore m4
margins, expression(predict(pr outcome(3)) + predict(pr outcome(4)))  at(lnpop=2.9957323 size_low=(75 150)) atmeans post pwcompare
mat list e(b) 


// Generate predicted effect of changes in party variables, as reported in paper:

foreach var of varlist hou_majority policy unified_gov {
quietly sum `var' if numb==1, detail
local rmin = r(min)
local rmax = r(max)
local rmean = r(mean)

estimates restore m2
margins, expression(predict(pr outcome(3)) + predict(pr outcome(4))) at ( `var'=`rmin' partyid=(-2 2)) atmeans post pwcompare
mat list e(b)

estimates restore m2
margins, expression(predict(pr outcome(3)) + predict(pr outcome(4))) at ( `var'=`rmean' partyid=(-2 2)) atmeans post pwcompare
mat list e(b)

estimates restore m2
margins, expression(predict(pr outcome(3)) + predict(pr outcome(4))) at ( `var'=0 partyid=(-2 2)) atmeans post pwcompare
mat list e(b)

estimates restore m2
margins, expression(predict(pr outcome(3)) + predict(pr outcome(4))) at ( `var'=`rmax' partyid=(-2 2)) atmeans post pwcompare
mat list e(b)
}

sum debt_net_pc_real
estimates restore m2 
margins, expression(predict(pr outcome(3)) + predict(pr outcome(4))) at(debt_net_pc_real = (0.0409679 6.453723)) atmeans
	

estimates restore m2
local est1 = ln(1.18)
local est2 = ln(.575)
margins, expression(predict(pr outcome(3)) + predict(pr outcome(4))) at(lncs = (`est1' `est2') )  atmeans

log close
